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Using both exact diagonalizations and diagonalizations in a subset of short-range valence bond 
singlets, we address the nature of the groundstate of the Heisenberg spin- 1/2 antiferromagnet on 
the square lattice with competing next-nearest and next-next-nearest neighbor antiferromagnetic 
couplings {J1 — J2 — J3 model). A detailed comparison of the two approaches reveals a region along 
the line (J2 + Jz)/ Ji — 1/2, where the description in terms of nearest-neighbor singlet coverings is 
excellent, therefore providing evidence for a magnetically disordered region. Furthermore a careful 
analysis of dimer-dimer correlation functions, dimer structure factors and plaquette-plaquette cor- 
relation functions provides striking evidence for the presence of a plaquette valence bond crystal 
order in part of the magnetically disordered region. 

PACS numbers: 75.10.-b, 75.10.Jm, 75.40.Mg 



I. INTRODUCTION 

Frustration can drive the low-energy physics of bidi- 
mensional Heisenberg quantum antiferromagnets very far 
from conventional semi-classical Neel-like phases. In such 
a case, the breakdown of long range magnetic order in the 
ground state leads the system to reorganize in a typical 
quantum state where only local antiferromagnetic corre- 
lations are present, namely a superposition of short range 
valence bond (SRVB) states. In this regime, the system 
opens a gap to the magnetic excitations and the SU(2) 
symmetry of the Hamiltonian is restored. However, in- 
side this general frame, the nature of the SRVB ground 
state (GS) can be very different from one system to an- 
other. In the simplest scenario the spatial symmetry of 
the Hamiltonian can still be broken, leading to a valence 
bond crystal phase (VBC) characterized by long range 
order in the dimer-dimer correlation function^. Alter- 
natively, all symmetries can be restored in a flat super- 
position of SRVB states to form a so-called spin hquid 
(SL). 

Far from being purely academic, the precise determina- 
tion of the GS nature is a crucial question in the context 
of quantum phase transitions^. For example, the "dccon- 
fined critical point" (DCP) scenario has been proposed 
as a new class of criticality to describe the Neel to VBC 
transition2i4. More importantly, the nature of the mag- 
netic background dramatically affects the holon/spinon 
(de)confinement properties of the corresponding doped 
systems. It is therefore believed to be a key ingredient 
to understand exotic metallic states. 

In practice, it is often hard to fully characterize the 
type, from crystal to liquid, of a SRVB phase. In this 
respect, one of the most archetypal example of such a sit- 
uation is the J1 — J2 Heisenberg 5 = 1/2 antiferromagnet 
on the square lattice, where frustration is controlled by 
the next nearest neighbor interaction J2. Despite many 



years of numerical and analytical efforts, no definitive 
picture emerged around the maximally frustrated point 
J2/J1 ~ 0.5, where the magnetic order disappears. The 
main point of this article is to introduce a general frame- 
work to study this kind of highly frustrated antiferro- 
magnet (the SRVB method) and to revisit the question 
on this specific model within an extended version of the 
Hamiltonian including a third neighbor J3 interaction : 
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At a classical \eve]~'^, the effect of frustration and 
competition between J2/J1 and J^j J\ leads to four or- 
dered phases described in figure ^ The effect of quan- 
tum fluctuations on this classical phase diagram is still 
an open question. In the past 15 years, the situation 
has been somewhat clarified for the pure J\ — J2 model, 
especially in a range of parameters far from the max- 
imally frustrated point J2/J1 ^ 0.5. For J3 = and 
J2/J1 ^ 0.4, the classical (tt, tt) Neel behavior is essen- 
tially conserved^ up to a small reduction of the stag- 
gered magnetization. On the other hand, for J2/J1 ^ 0.6 
an order by disorder mcchanismiS selects two collinear 
states at q = (tt, 0) and (0,7r). In the parameter range 
where frustration is the largest, 0.4 < J2/J1 < 0.6, the 
situation is much more involved. Beside the fact that 
many approaches (including spin- wave theory^, exact di- 
agonalizationsS, series-expansionii and large- A'^ expan- 
sionsi^) have now firmly established that quantum fluctu- 
ations destabilize the classical ordered ground state and 
lead to a quantum disordered singlet ground state with a 
gap to the flrst magnetic excitation, its precise nature is 
still controversial : a columnar valence bond crystal with 
both translational and rotational broken symmetriesi, a 
plaquette state with no broken rotational symmetrjji^ or 
even a spin-liquid with no broken symmctryi4 have been 
proposed (see figure QJ. 
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FIG. 1: (Color online) Classical phase diagram of the 
Ji ~ J2 — Js model : I Neel [11,11), II Collinear (0, tt) and 
(7r,0), III Helicoidal {q,q), IV Helicoidal {q,n). The snap- 
shots refer to the quantum version of the Ji — J2 model and in 
particular the various possible scenarios in the magnetically 
disordered gapped (red dotted line) phase around J2/J1 ~ 0.5 
: columnar VBC, plaquette VBC or spin liquid. VBC corre- 
lations are investigated in details along the black dashed lines 
in the present paper. 



For the J3 7^ case, as remarked by Ferreu^, the end 
point of the classical critical line (J2 + 2J^)/Ji = 1/2 
on the J3 axis is substantially shifted to larger values of 
J3 when quantum fluctuations are switched on. For the 
pure Ji — J3 model, in this region of large frustration, a 
non-classical (but still controversial) phase appears be- 
tween the Neel (tt, tt) and the spiral {q, q) phases : a 
VBC columnar statei^, a spin-liquidi& or a succession of 
a VBC and Z2 spin-liquid phasesii have been proposed. 
The complete phase diagram of the Ji — J2 — J3, quantum 
antifcrromagnct is expected to be even richer. Indeed, 
preliminary calculations^ pointed towards an extended 
region with a quantum disordered state. 

In this paper we investigate the maximally frustrated 
region of this phase diagram (J2 -I- J3)/ Ji ~ 1/2 (dashed 
line in figure^ using both exact diagonalizations and a 
SRVB method which consists in diagonalizing the Hamil- 
tonian in a subset of singlets states that can be written 
in term of SRVB states. In the first section, we intro- 
duce in details the method as a natural tool to study 
magnetically disordered phases and discuss its advan- 
tages and limitations. In the second part, we show nu- 
merical evidences for an extended non-magnetic phase 



around (J2 + J3) / Ji ^ 1/2. In a third part we present 
calculations and finite size analysis of dimer-dimer corre- 
lation functions and dimer structure factors that estab- 
lishes the existence of a s-wave plaquette ordered phase 
breaking only translational symmetry when J3 > J2 and 
(J2 + J3)/Ji ~ 1/2. This point is directly confirmed in 
the last part by an inspection of plaquettc-plaquctte cor- 
relations. We conclude by emphasizing the interest of 
the Ji — J3 model as an example of Neel to VBC quan- 
tum phase transition and discuss the implications of our 
results for the much debated Ji — J2 model. 



II. SRVB METHOD 

From a numerical point of view, investigating the low 
energy physics of 2d frustrated quantum antiferromag- 
nets is a difficult problem. Among the three well known 
high precision and controlled methods, two of them can- 
not be applied, at least for the moment : Density Ma- 
trix Renormalization Group (DMRG) is only efficient in 
one dimension and Quantum Monte Carlo (QMC) suffers 
from a severe sign problem on these systems. The third 
method, namely exact diagonalizations (ED), consists in 
a complete enumeration of the Hilbert space followed by 
an iterative solving of the eigenproblem. The main ad- 
vantages of this approach are (i) it is numerically exact, 
(ii) any observable is accessible (iii) spatial symmetries 
can be fully taken into account thus providing momen- 
tum resolved results. Unfortunately, the first step of the 
method faces the exponential growth of the Hilbert space 
with system size for finite available computing resources. 
Nevertheless, this method is still widely and successfully 
used and is the source of many firmly established results. 

However, if one compares highly frustrated quantum 
antiferromagnets to more conventional unfrustrated ones 
(typically Neel like) , a phcnomcnological review of known 
results shows that 

1. the role of the singlet sector is overwhelming at low 
energy due to the opening of a singlet-triplet gap, 

2. the breakdown of antiferromagnetic long range or- 
der favors the emergence of local singlet patterns. 

In this respect, it is tempting to build a more specific 
approach taking into account these two points in order 
to systematically reduce the Hilbert space to a relevant 
subset adapted to describe magnetically disordered sin- 
glet states. 

Following point 1, a first systematic reduction of the 
Hilbert space could be obtained by directly working in 
the singlet sector S ~ 0. Unfortunately, ED are not 
adapted to an explicit implementation of the SU(2) sym- 
metry of the Heisenberg Hamiltonian because, from a 
numerical point of view, eigenvectors of the total spin 
turn to be very complex objects for large systems. In 
practice, the Hilbert space used in ED is a set of 5^ = 
eigenvectors. The expected benefit of such a reduction 
would boiSi of order ^ 1/N. 
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(a) (b) 

FIG. 2: (a) Arbitrary range VB state, (b) Nearest neighbor 
VB state (NNVB). The oriented bond between two sites i and 
j stands for [i,j] = {1/V2){\ TUj) - | Uh))- 

In fact, a natural framework for fixed states has 
been developed years ago2£. Indeed, the whole singlet 
subspace can be generated using arbitrary range cover- 
ings of the lattice with VB states (see figure 121(a)) : 

i^>=n[*'ji' (2) 

where = (1/V2)(| Uij) - I Uh))- However, the 
practical relevance of these states is very limited because 
the number of dimer coverings for the complete graph 
is iV!/(2^/2(A^/2)!) - (iV/e)^/2 ^^ich is much larger 
than the size of the singlet subspace. As a direct con- 
sequence, this family of states is overcomplcte. Further- 
more it is certainly not specifically adapted to the de- 
scription of non-magnetic (quantum disordered) phases 
since any kind of singlet state, including the finite-size 
Neel state, could be constructed by an appropriate linear 
combination of arbitrary range VB states. 

Let us now examine point 2. A simple way to re- 
duce the number of coverings while keeping only short 
range correlations is to restrict the range of the dimers to 
short-range, for example nearest neighbor valence bond 
states (NNVB) (see figure |21 (b)). A general solution 
to the question of enumerating these states has been 
given by Fisher—. It is exponential ka^ for large N, 
with a w 1.34 (square lattice), a w 1.53 (triangular 
lattice) and a ~ 1.26 (kagome lattice). As expected, 
these numbers are much smaller than the total number 
of singlets thus providing the desired selection inside the 
singlet sector. Nevertheless, two important questions de- 
serve attention : (i) Are these states linearly independent 
? (ii) Which class of singlet states can be obtained by 
linear combinations of NNVB states ? The first ques- 
tion has not been addressed analytically but numerical 
calculations^S, show that, unless for very small systems 
on the triangular lattice, these states are linearly inde- 
pendent for the square, triangular and kagome lattices. 
Concerning the second question, it is clear that any state 
involving only short range spin-spin correlations, from 
VBC to SL, can be captured by SRVB states. On the 
contrary, Liang et al. showedS^ that magnetic long range 



order cannot be obtained from linear combinations of 
such configurations. 

As a partial conclusion, selecting a subset of SRVB 
states in the singlet space provides a convenient frame- 
work to study the low-energy singlet sector of highly frus- 
trated antiferromagnets. If the physics of a given prob- 
lem can be captured in this restricted basis, this kind 
of approach not only makes larger systems accessible to 
computation but also gives some insights about the na- 
ture of the GS ruling out any magnetic long range order. 
For technical details and illustrations of the method the 
reader can refer to previous publications^^iSSiSl. Never- 
theless, let us recall one of the most salient characteristic 
of the calculation : one crucial property of SRVB states 
is their non-orthogonality (see appendix). At a numeri- 
cal level, the problem of diagonalizing the Hamiltonian is 
then shifted to the so called generalized eigenvalue prob- 
lem (GEP): 

det {n - EO) = 0, (3) 

where O denotes the overlap matrix. The GEP, espe- 
cially when 7i and O are non-sparse matrices, cannot 
be efficiently solved iteratively. A rather time consum- 
ing complete diagonalization has to be performed, which 
makes use of spatial symmetries necessary for large clus- 
ters. 

Finally, let us remark that the GS computed with this 
method can be seen as the best variational approxima- 
tion of the exact GS using the restricted NNVB subset of 
states. However, even for a magnetically disordered exact 
ground state, the wave function almost certainly involves 
still finite range but more than only nearest neighbor VB 
states. As a consequence, this approach is not designed 
to provide the state-of-the-art variational approximation 
of the exact GS, but rather to capture in a small subset of 
physically suggestive states, the main part of the absolute 
ground state wave function neglecting finite range corre- 
lations refinements whose sole effect would be to slightly 
renormalize energies. In this respect, solving |j2J) is ex- 
actly equivalent to diagonalize a sophisticated effective 
Hamiltonian, namely the exact projection of the Heiscn- 
berg Hamiltonian on the chosen SRVB subspace. 



III. SRVB REGION 

One of the main drawbacks of the SRVB method is 
its lack of built-in control : solving 101 is always possi- 
ble even if the selected SRVB subspace is irrelevant to 
describe the low-energy sector of Ti. It is therefore nec- 
essary to make systematic comparisons between SRVB 
results and exact ones. 

To do so, let us consider an intermediate size cluster, 
namely A'' = 32, and compute the GS energy both by 
ED and NNVB diagonalizations, respectively -E™ and 
^NNVB r^Yie accuracy and thus the validity of NNVB 
approach can be tested by a measurement of the param- 



4 



O.7»""''"#-""10' 




0.6 • 



3.46% ^3.93% —4.59% ^5.55% ^^6.92%^K.80^^^B|.04S^^K85% 



J. 



^5.40% ^^6.77%^^ .2 



„ „ 1.72% 2.18% 2.70% 3.37% ^ 4.25% .^5.40% ^6.77%^^.25% 

0.5 • • • • « ^ 



„ , 0.73% 0.88% 1.42% 2.18% 3.12% ^ 4.32% ^5.83% ^fc7.61% 

0.4 • * * 9 V 



0.3 



0.2 



0.1 



1.09% 0.51% 0.51% 1.18% 2.32% ^3.79% ^ 5.56% .^7. 



2.32% 3.79% ^5.56% .^7.1 

ITT? 



2.44% 1.32% 0.59% 0.52% 1.37% 3.13% ^5.4 



^4.32% 2.93% 1.73% 0.90% 0.83% 2.02% _ 4.82%^^fc.01% 



.6.44% _4.97% ^3.55% 2.28% 1.42% 1.48% 3, 



0.1 0.2 0.3 0.4 0.5 0.6 



9I^^|^^7' 



FIG. 3: Systematic comparison of ED and NNVB ground 
state energy for A'^ = 32. The radius of the circles is pro- 
portional to the NNVB ground state accuracy (Eq^^^ — 
Eo°)/Eo°. Typical values of the energies are given in ta- 
ble |T| 



eter (£;NNVB _ ^eD)/^ed_ q^^ figureEl this quantity is 
plotted as a function of J2/J1 and J3/J1. 

As expected, the NNVB ground state fails to approxi- 
mate the exact one in the regions of the phase diagrams 
known to be magnetically ordered : ( J2 <C Ji, J3 ^ Ji), 
( J2 » Ji, J3 < ^1) or (J3 » Ji, J2 < Ji). On the oppo- 
site, in the highly frustrated regime, an extended region 
of the phase diagram emerges around (J2 + Jz)/ Ji '--^ 1/2 
where (£;NNVB _ ^ED)/^ED gj^allcr than 1.5% and as 
small as 0.5% (see figure|21and tableP). 



J3) 




gNNVB 


(J2, J3) 


BED 


^NNVB 


(0.0,0.3) 


-18.71704 


-18.51215 


(0.2,0.4) 


-16.66878 


-16.43163 


(0.0,0.4) 


-18.12399 


-17.99224 


(0.3,0.1) 


-17.12863 


-16.97424 


(0.0,0.5) 


-17.92509 


-17.61700 


(0.3,0.2) 


-16.50461 


-16.41946 


(0.1,0.2) 


-18.43435 


-18.19099 


(0.3,0.3) 


-16.17630 


-15.98600 


(0.1,0.3) 


-17.72089 


-17.63119 


(0.4,0.0) 


-16.90813 


-16.66731 


(0.1,0.4) 


-17.33094 


-17.17892 


(0.4,0.1) 


-16.21783 


-16.08331 


(0.2,0.2) 


-17.39604 


-17.29400 


(0.4,0.2) 


-15.80152 


-15.58522 


(0.2,0.3) 


-16.86835 


-16.78183 


(0.5,0.0) 


-16.00307 


-15.76633 



TABLE I: ED and NNVB ground state energy for = 32 as 
a function of (J2, J3) (units of Ji). 

Before going any further in the analysis, it is important 
to have in mind the order of magnitude of the NNVB 
truncation of the Hilbert space. For such a system size, 
the dimension of the GS representation (k = (0,0), s- 
wave) is 1184480. This has to be compared to the number 
of NNVB configurations in the same representation which 



is only 182. The reduction factor is thus ^ 10^. 

Considering both the accuracy of Eq^^^ and the 
rather drastic reduction of the singlet space we can con- 
clude to the existence of an extended region in the phase 
diagram, around {J2 + J3)/Ji ^ 1/2, where the exact 
GS can be described with only NNVB states. Never- 
theless, in order to investigate the precise nature of the 
ground state using this wave function, it is important 
to go beyond this energetic criterion. A direct evalua- 
tion of the overlap between the exact ground state and 
the NNVB variational wave function {ipoli^o^^^} can- 
not be done easily, but it is straightforward to com- 
pute an upper bound for the so-called "missing weight" 
1 — |(V'o|V'o"^^^)P which, crudely, quantifies the "accu- 
racy" of the wavcfunction w.r.t. the exact GS. A formal 
normalized expansion of IV'^^^^) = X^i'^ilV'i) on the 
exact eigenstates leads to the expression of Eq^^^ = 
la^pi^i as a function of the exact eigenenergies Ei. 
Since Ei > Ei for z > 1 one obtains, 
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(4) 



This quantity is represented on figure ^ as a function of 
J2/J1 and J3/J1. Despite the fact that this upper bound 
is far from being optimal since Ei is only a crude lower 
bond for highly excited states, the same region of the 
phase diagram (as the one determined previously on a 
purely energetic criterion) emerges where KV'olV'o"^^^)! 
is at least 90% in the worst case and up to 95% in the 
best case. 

This picture clearly confirms that around (J2 -I- 
•h)/Ji ~ 1/2 the essential part of the GS wave func- 
tion can be captured using only few SRVB states, namely 
NNVB configurations. As mentioned in the previous 
section, there is no doubt that this accuracy could be 
systematically improved by dressing IV'?'^^^) with some 
longer (but still finite) range VB configurations (eg next 
nearest neighbor VB configurations). Although includ- 
ing such additional configurations are expected to lower 
even further the variational energy, this would be no more 
than refinements and we believe that the approach here 
already fully captures the physical picture of a SRVB 
ground state. 



IV. DIMER-DIMER CORRELATIONS AND 
STRUCTURE FACTORS 

The next important question is now to investigate the 
nature, VBC or SL, of the SRVB ground state in this re- 
gion. To address this question we used the SRVB method 
to compute the dimer-dimer correlation function : 



C, 



ijkl 



((S,.S,)(S,.S0)-((S,.S,))2 



(5) 



The SRVB method allows a systematic computation of 
(0) on the extended SRVB phase for cluster sizes ranging 
from TV = 20 to TV = 40. 




FIG. 4: (Color online) Upper bound for the "missing weight" 
1 - K^o|V'S"~^^^)r The radius of the blue circles is pro- 
portional to the upper bound given in equation Q. Values 
greater than one being irrelevant are represented as unit ra- 
dius red circles. 



Real space picture. Figures [S] and |^ arc snapshots of 
the resuhs for = 40 respectively for the pure Ji — J2 
model at J2/J1 = 1/2 and the pure Ji — J3 model at 
J3/J1 = 1/2. Both systems exhibit, for bonds parallel 
to the reference bond a clear alternating pattern 

of correlated and anticorrclatcd rows. Moreover, around 
the maximal distance from the reference bond, the values 
of the parallel correlations are almost constant. As a 
consequence both figures [S] and suggest a translational 
symmetry breaking VBC phase with a stronger signal in 
the latter case. 

As suggestive as this kind of picture may be, two im- 
portant question have to be addressed : (i) To what kind 
of VBC phase figures [S] and correspond ? (ii) Is this 
suggested long range order robust when N 00 7 

Even if at first sight these real space pictures naively 
suggest a columnar arrangement, a plaquette VBC order 
cannot be ruled out. In order to investigate the nature 
of the VBC ground state, we introduce 3 trial wave func- 
tions ipci "fps and ipd respectively referring to a colum- 
nar, s-wave plaquette and d-wave plaquette state (see 
Appendix). These wave functions are designed to have 
the same symmetry as the finite size ground state, namely 
k = (0, 0) s-wave, in order to allow direct comparisons 
with the numerical results. 

The computation of the dimer-dimer correlations in 
these wavefunctions is presented in details in the Ap- 
pendix and the results are summarized in table llTTl First, 
a comparison of our previous numerical results with those 
of table mil shows that the d-wave plaquette scenario is 
very unlikely. Furthermore, the results of the Appendix 
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FIG. 5: (Color online) Dimer-dimer correlation function 
for a 40 site cluster with periodic boundary conditions with 
J2/J1 = 1/2 and J3 = 0. The dashed line delimits the clus- 
ter, (i,j) is the reference bond, and the width of the solid 
bonds {k,l) are proportional to the absolute values of Ciju- 
The blue (resp. red) bonds denote positive (resp. negative) 
correlations. Numbers correspond to |10''Cijfe;j. 
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FIG. 6: (Color online) Same as figureOfor J3/J1 = 1/2 and 
J2 =0. 



suggest that the key criterion to discriminate between a 
pure columnar and a pure s-wave plaquette VBC, on the 
basis of dimer-dimer correlations, is the ratio between (i) 
perpendicular bond correlations (with respect to the ref- 
erence bond) and (ii) parallel bond correlations in odd 
columns (defining the reference bond column as even). 
In the first case it is expected to be equal to 1 while it 
should vanish in the latter case. 

For the data shown in figures |S1 and |3 if one considers 
the most distant bonds from the reference one, the typical 
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FIG. 7: (Color online) Phase factors e\{k,l) for structure 
factors : (a) "VBC" (b) "col". Dashed bonds are either in- 
cluded or excluded form the definition @ in the fitting pro- 
cedure in order to test the sensitivity of the extrapolation 
scheme to irrelevant short range contributions. Note that the 
(fc, I) bonds nearest neighbors to the reference one (central 
black solid bond) are always omitted in the sum defining Sx. 



value of this ratio is of order 1/20 and 1/100 respectively. 
This strongly supports a s-wave plaquettc scenario for 
J3 = Ji/2 and J2 = while the situation appears more 
involved for J2 = Ji/2 and J3 = where the ratio is 
still very small but for a much weaker overall long range 
correlation signal. 

Finite size analysis and structure factors. It is crucial 
to study the robustness of this picture with the system 
size. A convenient way to investigate the thermodynamic 
limit is to introduce spatially integrated quantities such 
as dimer structure factors and perform finite size scaling. 
The essential difference between columnar and s-wave 
plaquette orders is the breakdown of rotational symme- 
try. Following Ref. it is possible to build two structure 
factors 5'vBC smd Scoi with the following properties : 

• SyBC diverges at thermodynamic limit both in 
columnar and plaquette states, 

• 5coi diverges at thermodynamic limit only in a 
columnar state. 

To achieve this, the form factors introduced in S'vbc a-nd 
5coi have to reflect the patterns of table IIIll It is easy to 
verify that appropriate structure factors can be defined 
e.g. as : 



S\ = y^^e\{kJ)Cijki, 



(6) 



where A stands for cither "VBC" or "col" and the cor- 
responding form factors e\{k,l) are defined according to 
figure [7| 

In an ordered phase 5a is extensive so that Sx/Ni,, 
where Nb denotes the number of bonds involved in lO, 
is expected to scale like + A/N with being the 
square of the bond order parameter in the thermodynamic 
limit. The divergence (resp. finite value) of Sx is thus 
signaled by a finite (resp. vanishing) . 
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FIG. 8: (Color online) Left panel (a) : SyBc/Nt as a function 
of l/N along the J2 + J3 = Ji/2 line. Note that the N = 36 
data is excluded from the linear fits represented as dashed 
lines. Right panel (b) : Same as left panel with a modified 
definition of S'vbc in which very short range contributions are 
excluded (dashed bonds in figure |7|l. 



We performed this type of scaling on S'vbc /N^b for N = 
20,32,36 and 40 along the line (J2 + Ja)/^! = 1/2. As 
shown in figurc|Hl(a), the quality of a 1/A'' extrapolation 
is greatly affected by the N = 36 data. This point is 
due to the peculiar shape of this cluster whose periodic 
boundary conditions induce short loops that have the 
tendency to overestimate the influence of the reference 
bond and thus Sx . We therefore excluded this set of data 
in the analysis depicted on figure |Sl(a). Along the whole 
(J2 + •h)/Ji ~ 1/2 hue, the fit reveals a non- vanishing 
extrapolated and a standard evaluation of errors 

bars on the extrapolated values is presented on figure |H1 
(b) (thin line labelled "No cut"). 

/^From a technical point of view it is fair to evaluate, in 
the extrapolation scheme, the influence of the strong con- 
tributions to the structure factor coming from the short 
range part of the dimer-dimer correlations (see figures 
and El . There are at least two reasons to discuss this 
aspect : (i) The short range part of the data is irrel- 
evant at large distance and therefore, a non negligible 
contribution to the thermodynamic extrapolation would 
indeed be problematic (ii) As shown in the Appendix, a 
substantial enhancement of the short range dimer-dimer 
correlations is expected to occur in plaquette states (see 
table Im)) . 

The sensitivity of the fit to the (irrelevant) short range 
correlations can be tested by systematically removing 
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FIG. 9: (Color online) Left panel (a) : Extrapolation of 
SvBc/Nb as a function of J2/J1 along the J2 + J3 = Ji/2 line 
(Thick solid line with error bars). Finite size Svbc/Ni, data 
for A'^ = 20, 32 and 40 are represented as thin lines and circles. 
The error bars reflect the quality of the 1/A^ fit presented on 
figure|S](b). Right panel (b) : Influence of short range contri- 
butions to the extrapolated VBC structure factor Cybc along 
the line (J3 + J2) = 1/2 as a function of J2/J1 and compar- 
ison with expected values for columnar and s-wave plaquette 
states. Thin (resp. thick) line with error bars labelled "Not 
cut" (resp. "Short range cut") corresponds to the results of 
the fits of SvBc/Nb including all range contributions (resp. 
excluding short range contributions) represented on figure |H| 
(a) (resp. (b)). Thick dashed lines are the expectations val- 
ues of the structure factors C^^c at thermodynamic limit for 
the pure columnar (short dashed line) state ipc and the pure 
s-wave plaquette state tps- 



from the sum defining the structure factor the contri- 
bution of the neighboring bonds of the reference one (see 
dashed bonds in figure [Tj). As shown in figures |H1 and |H| 
(b), when J2/J1 1/2 the extrapolated values of S'vbc 
is insensitive to the short range correlations, while in the 
crystalline phase, the procedure of removing the short- 
range part of the data has a systematic tendency to en- 
hance the VBC order parameter and to lower the error 
bars thus improving the confidence of the extrapolated 
value. This fact convincingly establishes that the under- 
lying GS has a VBC long order for J2/J1 < 0.2 - 0.3 
but also gives some further indication: very short range 
dimer-dimer correlations in the GS are responsible for a 
slight perturbation of the extrapolation which is compat- 
ible with the local enhancement of Cijki observed in the 
trial plaquette state ^pg when {k,l) is lying next to 
(see table mH) . From a technical point of view, in order 
to exclude this kind of short range effect, we exclude for 
further analysis the short distance contribution to the 
definition ((HJ of S'a. 

A careful inspection of figure |H1 reveals two regimes of 
parameters for J2/J1 : below ~ 0.2 the opening of the 
errors bar is due to to a convex deviation from a perfect 
linear behavior, while it is concave above J2/J1 ~ 0.3. 
As a consequence, the extrapolation scheme respectively 



underestimates and overestimates S'vbc- This confirms 
that the crystalline order is indeed robust for J2/J1 < 
0.2 — 0.3. Moreover the extrapolated value for J2/J1 = 
(and J3/J1 = 1/2) is 0.032 ± 0.003 which compares very 
well with the expected values of the pure columnar or 
plaquette crystalline states which respectively equal to 
9/256 - 0.035 and 1/32 - 0.031 (see table |llll in the 
Appendix and figure |Sl(b)). 

In contrast, due to large error bars and the slight con- 
cavity of the 1/A'' extrapolation, a vanishing Cysc 
not be ruled out from our data for J2/J1 larger than 0.3 
and therefore the existence of a crystalline long range 
order for J3 = is not proven by the present calculation. 

Let us now turn to S'coi. The size dependence of 
Scoi/Nh does not allow a confident extrapolation to ob- 
tain C^j with enough accuracy. Nevertheless, for all clus- 
ters S'coi is always a very small fraction of S'vbc SiS can 
be seen by comparing figures (a) and ^] for iV = 32 
and N = 40. Typically the ratio Scoi/Svbc is of or- 
der 1/20 for J2/J1 = and 1/15 for J3/J1 = 0. The 
expected values of this ratio for the pure columnar and 
s-wave plaquette state (see table IIII|) are respectively 1 
and 0. 

We cannot draw definitive conclusions from our data 
in the regime where J2/J1 ^1/2 and J3 ^ since our 
scaling does not exclude a scenario where CyB^ and 
would vanish. In contrast, on the {J3-\- J2) / Ji = '^l'2, line 
for small J2 and up to J2/J1 ^ 0.3, the fact that S'coi 
is much weaker than S'vbc is very much in favor of the 
s-wave plaquette scenario with an absence of rotational 
symmetry breaking and seems to rule out a simple long 
range columnar order for which Scd ~ S'vbc in the ther- 
modynamic limit. Note that a small spatial anisotropy of 
the plaquette phase is still possible. This scenario where 
the vertical and horizontal bond amplitudes within the 
resonating plaquettes are slightly different would indeed 
lead to a small value of the columnar structure factor 
in the thermodynamic limit and a GS degeneracy of 8 
(instead of 4)2^. 

V. PLAQUETTE-PLAQUETTE 
CORRELATIONS 

A careful analysis of the difference in dimer-dimer cor- 
relations in a columnar dimer versus an s-wave plaquette 
ordered singlet state performed in the previous section 
yielded strong support for a plaquette phase. In order to 
directly image the plaquettes in real space we calculate 
the following 8-spin correlation function using ED: 

Cpiaquettesb, g) = {QpQq) - {Qp)"^ (7) 

where p and q denote two different plaquettes and ^ 
denotes the cyclic exchange operator of the four spins 
on a given plaquette. This correlation function has also 
been used in a recent study of plaquette order in the 
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FIG. 10: Comparison between SyBc/Ni, and Scoi/A^b as a 
function of J2/J1 along the J2 + J3 = Ji/2 line. Thin lines 
with circles : finite size data for Xcoi/Nt for A'' = 32 and 
N = 40. Thick dashed lines are the expectations values of 
the structure factors Cvbc C^i at thermodynamic limit 
for the pure columnar (short dashed line) state tpc and the 
pure s-wave plaquette state Vs- Note that C^i = in the 
s-wave plaquette state. The thick line with errors bars is the 
same as in figure |3] 



checkerboard antiferromagnet2&. If we want to discrim- 
inate between a columnar dimer state and a plaquette 
ordered state in the following, it is useful to note that 
in a columnar dimer state one has two distinct expecta- 
tion values of {Qp) (either covering two singlet bonds or 
none), whereas in a plaquette ordered state we expect 
three distinct expectation values (on a singlet plaquette, 
between two adjacent singlet plaquettes, or sharing the 
corners of four distinct singlet plaquettes). This num- 
ber is expected to translate into the number of different 
values in the correlation function Eqn. |(7J). 

We present the results obtained by ED on a = 32 
sample in Fig. ^] both along a line with {J2 + 'h) / Ji = 
1/2 (upper row) and along the pure J3 line (lower row). 
In the cases where strong correlations are seen, we basi- 
cally detect three different types correlation function val- 
ues, in agreement with the expectations of the plaquette 
phase, as pointed out above. Furthermore the spatial 
structure coincides with the plaquette picture, i.e both 
the positively and the negatively correlated plaquettes 
form a distinct 2x2 superlattice, shifted by the vec- 
tor (1,1) with respect to each other. The evolution of 
the correlations as a function of J2 and J3 shows that 
the strength of the correlations both decreases as one 
moves away from the point J2 = 0, J3/J1 = 1/2 ei- 
ther along the pure J3 line or along the line with fixed 
(J2 + 'h)/ Ji = 1/2, in agreement with the results of the 
preceding section based on dimer-dimer correlations. In- 
terestingly the correlations at the much debated point 
J2I J\ = 1/2, J3 ~ are rather weak, but still carry some 
remnants of the plaquette phase, at least for this iV = 32 
sample. 



VI. DISCUSSION AND CONCLUSIONS 

An extensive numerical study of the Heiscnberg Ji — 
^2~</3 antiferromagnet using both exact diagonalizations 
and a short range valence bond method shows that, in 
the most frustrated part of the phase diagram (around 
J2 + J3 ~ t^i/S), the ground state can be captured us- 
ing only nearest neighbor valence bond coverings of the 
square lattice. The emergence at low energy of short 
range valence bond singlet physics for these parameters 
and thus the breakdown of magnetic long range order is a 
direct consequence of the strong frustration of the model. 
Moreover, we characterize the ground state by an anal- 
ysis of dimer-dimer correlations, dimer structure factors 
and plaquette-plaquette correlations and show numeri- 
cal evidences for an extended valence bond crystal phase 
around J2 + J3 — Ji/2 and J2 < J3 where the ground 
state is a s-wave plaquette state only breaking transla- 
tional symmetry. As a consequence, the Ji — J3 model 
provides an example of frustration-driven Neel to VBC 
quantum phase transition. Note that the SRVB frame- 
work can be readily extended to include singlet pairs be- 
yond nearest neighbors. However, we believe that this 
will modify only slightly the results in the maximally 
frustrated region where the magnetic correlation length 
is very small. Such an approach could nevertheless be 
useful to investigate properties close to the critical point 
where the spin correlation length is expected to grow. In 
that respect such a transition can be probed by intro- 
ducing static (non magnetic) impurities^. Again, our 
framework could be extended to that case^i. 

For J3 < J2, including the much debated frustrated 
phase of the pure Ji — J2 model, the NNVB description 
of the ground state remains relatively robust. While our 
results are not able to resolve the controversy around 
J2/J1 « 1/2, the inclusion of an additional J3 coupling 
allows us to put this region into a broader perspective. 
Wc show that an antiferromagnetic J3 is useful in push- 
ing the magnetically ordered phases further apart, there- 
fore leaving more room for the disordered phases, and 
enabling us to reveal a robust plaquette singlet ordered 
phase. On the contrary, a ferromagnetic J3 interaction 
will probably lead to a direct first order transition be- 
tween the (tt, it) and the (tt, 0) Neel order phases as func- 
tion of J2 , similar to the classical analysis and numerical 
results on the related bcc lattice^S. The closeness of the 
magnetically order phases and the related phase transi- 
tions are probably responsible for the enormous difficulty 
in settling the controversy on the nature of the magnet- 
ically disordered phase (s) of the pure Ji — J2 model on 
the square lattice. 
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FIG. 11: Plaquette correlation function Cpiaqucttos(p, ?) [Eqn. (|7|l], obtained by exact diagonalization on a TV = 32 sample. 
The black squares denote the reference plaquette, filled blue circles correspond to positive values and empty red circles denote 
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APPENDIX: VB STATES PROPERTIES AND 
DIMER-DIMER CORRELATIONS FOR SOME 
RELEVANT VBC STATES 
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In this appendix, we recall some basic overlap proper- 
ties of VB states and compute dimcr-dimcr correlations 
expectation values for columnar and plaquette states. 

Overlaps. Two VB states \ip) and \ip) have an non 
vanishing overlap {(p\ip). To compute this quantity it 
is convenient to consider the loop diagram obtained by 
superimposing both configurations (see. figure [T^ . Be- 
cause loops are decoupled, {(p\^) is the product of each 
loop contribution. Since there is only two ways to de- 
scribe any loop with antiparallcl spins, | ti ... I) and 
I It ■•■ the overlap is 2"' (up to a normalization 
constant) with ni the total number of loops. The nor- 
malization is fixed by {ip\(p) = 1 which diagram contains 
N/2 trivial loops. The result is then {iplip) = £^,.^2"'-^/^ 
where the sign e^^^ is due to the relative orientations of 
dimers in \ip) and In the case of nearest neighbor 
VB on a bipartite lattice this sign can be fixed to 1 by 
convention, but in general this sign cannot be considered 
as constant. 

Orthogonal states in the thermodynamic limit. Let us 
consider two VB states such that the loop diagram con- 
tains at least one large loop, namely a loop involving 
aN^ sites with /? 7^ 0. The number of remaining sites 
is — aN^ so the maximal total number of loops is 
1 + {N -aN'^)/2. Hence, \{ip\^) \ < 2i-(i/2)«^'' and the 



FIG. 12: Overlaps between two VB states. 



overlap goes to zero when N goes to infinity. 

Another class of orthogonal states at thermodynamic 
limit is formed by states whose loop diagram contains an 
extensive number of loops ni = aN (note that a < 1/2). 
If a < 1/2, then |(<p|-0)| < 2^("-i/2) and the two states 
arc orthogonal when N goes to infinity. 

Columnar state and plaquette states. We define the 
columnar state (resp. plaquette state) as the equal 
weight linear combination of the four states (see figure 
I13II obtained by translation of the columnar (resp. pla- 
quette) covering of the lattice. The resulting state has a 
k = (0, 0) momentum, thus allowing direct comparisons 
with the finite size k = (0, 0) GS discussed in the article. 
Note that two different plaquette states can be defined 
on four sites : one is symmetric upon rotation of the pla- 
quette (see figure PHI a) and the second is antisymmetric 
(see figure ITHb'). Wc refer to these states respectively as 
■0s and ipd- The columnar state is denoted by -Re- 
using the results of the previous paragraph we can 
show that the four components of ^pc are mutually or- 
thogonal in the thermodynamic limit. It is easy to check 
(sec figure 1131 top row) that the overlap diagram have 
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FIG. 13: Definitions of columnar (top row) and plaquette 
(bottom row) state. 




■ • • a 

a 

FIG. 14: Definitions of s-wave plaquette (a), d-wave plaque- 
tte (b) and inequivalent couples of bonds {a, 13), (0,7), {ct,S) 
(c). 

either at least one large loop (in fact ~ \/N) or an ex- 
tensive number of loops. 

The very same argument can be applied for the four 
components of ipd after rewriting each d-plaquette as 
crossing dimers along the diagonals (see figure [HI bV 

The case of tps deserve more attention. Let us define an 



operator U that change the orientation of all the dimers 
on half of the vertical (or horizontal) lines in an alter- 
nating pattern. This operator is trivially self-adjoin and 
= id, so is a unitary transform and thus conserves 
scalar product : {(p\ip) = {fu\4'u) with \^pu) = l^W) ^-nd 
\ijju) ~ ^1^)- Since the action oiU on a, plaquette cov- 
ering simply exchange s-plaqucttes and d-plaquettes (see 
figures [T^ a and b), the orthogonality of the four compo- 
nents of is shown. 

The absence of interference between components of 
t/'cV's or V'd 9'lso occurs in the computation of {Pa) or 
{PaP^i) where Pb denotes the operator that permutes the 
two sites of bond h and ^ = (5,^,5 (see figure H"!) b). 
Indeed, the permutation of two or four sites on one com- 
ponent does not affect the existence of either ^ \fN large 
loops nor an extensive number of loops when overlapping 
with another component. 

Dimer-dimer correlations. The aim of this section is to 
compute {PaPfi) - {Pa)"^ ■ Note that since = 2Si.Sj + 
1/2, the correlation (PqP^) — {Pa)'^ is just related to the 
same expression with spin operators by a factor 4. 

By a direct evaluation we derive the basic rules to com- 
pute {Pa) for one component of ij^c, V's or tp^ : 

• {Pa)c = — 1 if a is occupied by a dimcr, +1/2 oth- 
erwise, 

• {Pa)s = ^1/2 if ct belongs to a plaquette, 4-1/2 
otherwise, 

• {Pa)d = +1/2 whatever a belongs or not to a pla- 
quette. 

Using these rules it is possible to evaluate {PaPfj.) — 
{Pa)'^ by a simple inspection of the 4 components contri- 
butions (see figure EJ) of i/jc, '(ps or ipd as shown in table 
UTI We summarize in table lllll the expected dimer-dimer 
correlation values in units of permutations and spin oper- 
ators as well as expected VBC and Columnar structure 
factors according to definition ©. Note that we only 
consider bonds /3, 7 and 5 that do not share sites with a. 
Also remark that for plaquette states {tps and ipd), very 
short range (a,/3) correlations differ from longer range 
ones when (a,/3) belong to the same plaquette. These 
short range anomalies are reported in tables ITTl and lllll in 
italic. 



1 N. Read and S. Sachdev, Phys. Rev. Lett. 62,1694 (1989). 

^ S. Sachdev, Quantum Phase Transitions (Cambridge Uni- 
versity Press, Cambridge, 1999). 

^ T. Senthil, A. Vishwanath, L. Balents, S. Sachdev and 
M. P. A. Fisher , Science 303, 1490 (2004). 

* T. Senthil, L. Balents, S. Sachdev, A. Vishwanath and 
M. P. A. Fisher, J. Phys. Soc. Jpn. 74 SuppL, 1 (2005). 

^ A. Moreo, E. Dagotto, Th. Jolicoeur and J. Riera, 
Phys. Rev. B 42, 6283 (1990). 



^ A. Chubukov, Phys. Rev. B 44, 392 (1991). 

J. Ferrer, Phys. Rev. B 47, 8769 (1993). 
^ P. Chandra and B. Dougot, Phys. Rev. B 38, R9335 (1988). 

E. Dagotto and A. Moreo Phys. Rev. Lett. 63, 2148 (1989); 

D. Poilblanc, E. Gagliano, S. Bacci, and E. Dagotto Phys. 

Rev. B 43, 10970 (1991); H.J. Schulz, T. Ziman and 

D. Poilblanc, J. Phys. I (France) 6, 675 (1996). 
" P. Chandra, P. Coleman and A.I. Larkin, Phys. Lett. 64, 

88 (1990). 



11 



Trial state 


Columnar (xpc) 


s-wave Plaquette 




d-wave Plaquette 




pairs of bonds 




(a, 7) 


{a,S) 


(a 




(",7) 


{a,S) 


(a 


/3) 


(q,7) 


{a,S) 


(a) 


+ 1 


-1/2 


-1/2 


+ 1/4 


+1 


-1/4 


+ 1/4 


+ 1/4 


+1 


+ 1/4 


+ i/4 


(b) 


+1/4 


-1/2 


+1/4 


+ 1/4 


+1/4 


-1/4 


-1/4 


+ 1/4 


+1/4 


+ 1/4 


+ 1/4 


(c) 


+1/4 


+1/4 


-1/2 


+ 1/4 


+1/4 


-1/4 


-1/4 


+ 1/4 


+1/4 


+ 1/4 


+ 1/4 


(d) 


+1/4 


+1/4 


+1/4 


+ 1/4 


+1/4 


-1/4 


+1/4 


+ 1/4 


+1/4 


+ 1/4 


+ 1/4 


Mean value 


7/16 


-1/8 


-1/8 


+ 1/4 


+7/16 


-1/4 





+ 1/4 


+7/16 


+ 1/4 


+1/4 



TABLE 11: (PaPb) — (Pa)^ with b = /3, 7,(5 computed in eacii of tiie four components (labelled from (a) to (d) in figure ITHll of 
the 3 trial states V'c, V'a xpd- Italic values correspond to the peculiar short range case where (a, 13) share the same plaquette. 



Trial state 


Columnar (^c) 


s-wave Plaquette (^s) 


d-wave Plaquette (V'd) 




1/8 





1/2 


{PcPp) - {Pcf 

((S.S).(S.S)^) -{(S.S)„)2 
Normalized (a, /3) 


7/16 
27/64 
27/256 
1 


1/4 
1/4 
1/16 
1 


7/16 
7/16 
7/64 
7/4 


1/4 





7/16 
3/16 
3/64 
3/7 


(Pc^P-y) 
{Pc^P-,} - {Pcf 

((S.S)a(S.S)^) -((S.S)„)2 
Normalized (a, 7) 


-1/8 
-9/64 
-9/256 
-1/3 






-1/4 
-1/4 
1/16 
-1 


1/4 





(PaPs) 
(P Px) - (P 

((S.S)..(S.S)0> -{(S.S)c.)2 
Normalized {a, 5) 


-1/8 
-9/64 
-9/256 
-1/3 








1/4 





Correlation snapshots 
































































































































































'^VBC 


9/256 = 0.035156 


1/32 


= 0.03125 





^col 


9/256 = 0.035156 








'-'col/'^VBC 


1 





Undefined 



TABLE IIL Expectations values of correlations and structure factors for V'c, tp3 and ^pd- Note that for the plaquette states, the 
correlations on the bonds next to reference one differs from the others (see italic numbers in columns 3 and 4). 



^1 M.P. Gelfand, R.R.P. Singh and D.A. Huse, Phys. Rev. B 
40, 10801 (1989); V. N Kotov, J. Oitmaa, O. Sushkov and 
Zheng Weihong, Phil. Mag. B 80, 1483 (2000); O. Sushkov, 
J. Oitmaa and Zheng Weihong, Phys. Rev B 66, 054401 
(2002). 

N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991). 
M.E. Zhitomirsky and K. Ueda, Phys. Rev. B 54, 9007 
(1996); L. Capriotti and S. Sorella, Phys. Rev. Lett. 84, 
3173 (2000); K. Takano, Y. Kito, Y. Ono and K. Sano, 
Phys. Rev. Lett. 91, 197202 (2003). 

L. Capriotti, F. Becca, A. Parola, and S. Sorella, Phys. 
Rev. Lett. 87, 097201 (2001). 

P.W. Leung and N.W. Lam, Phys. Rev. B 53, 2213 (1996). 
^® L. Capriotti, D.J. Scalapino and S.R. White, 
Phys. Rev. Lett. 93, 177004 (2004). 

L. Capriotti and S. Sachdev, Phys. Rev. Lett. 93, 257206 
(2004). 



^® F. Figueirido, A. Karlhede, S. Kivelson, S. Sondhi, M. Ro- 
cek, and D. S. Rokhsar Phys. Rev. B 41, 4619-4632 (1990). 
For N spin 1/2, the sizes of = and 5*" = sub- 
spaces are : A/'gs^o = N\/{{N/2y.{l + N/2)\) and A^s-^o = 
N\/iiN/2y.{N/2y.). 

^° L. Hulthen, Ark. Mat. Astron. Fys. A 26, (11), 1 (1938); 
M. Karbach, K.-H. Miitter, P. Ueberholz and H. Kroger, 
Phys. Rev. B 48, 13666 (1993). 

21 M. E. Fisher , Phys. Rev. 124, 1664 (1961). 
M. Mambrini, unpublished. 

S. Liang, B. Dougot and P. W. Anderson, Phys. Rev. Lett. 
61, 365 (1988). 

'^^ M. Mambrini and F. Mila, Eur. Phys. J. B 17, 651 (2000). 
On a finite cluster, this degeneracy is removed and leads 
to quasi-degenerate eigenstates of different quantum num- 
bers. 

J.-B. Fouet, M. Mambrini, P. Sindzingre, and C. LhuiL 



12 



lier, Phys. Rev. B 67, 054411 (2003); see also S.E. Palmer 
and J.T. Chalker, Phys. Rev. B 64, 94412 (2001) and 
B. Canals, Phys. Rev B 65 184408, (2002). 
S. Dommange, M. Mambrini, B. Normand and F. Mila, 
Phys. Rev. B 68, 224416 (2003). 



D. Poilblanc, A. Lauchh, M. Mambrini, and F. Mila Phys. 
Rev. B 73, 100403 (2006). 

R. Schmidt, J. Schulenburg, J. Richter, and D. D. Betts, 
Phys. Rev. B 66, 224406 (2002). 



